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Abstract. The techniques which allow the numerical evaluation of the 
statistical properties of the potential energy landscape for models of simple liquids 
are reviewed and critically discussed. Expressions for the liquid free energy and 
its vibrational and configurational components are reported. Finally, a possible 
model for the statistical properties of the landscape, which appears to describe 
correctly fragile liquids in the region where equilibrium simulations are feasible, 
is discussed. 



PACS numbers: PACS number(s): 

1. Introduction 

Understanding the dynamic and thermodynamic properties of supercooled liquids is 
one of the more challenging tasks of condensed matter physics ( jj], H> H @> H> H ) • ^ 
significant amount of experimental [0, ||, ||, |l0| , numerical , and theoretical work 
|{l2[ |l3| , p"4[ Hp , is being currently devoted to the understanding of the physics of 
the glass transition and to the associated slowing down of the dynamics. In the latest 
years the study of the topological structure of the potential landscape (PEL) [jlTj and 
the connection between the properties of the PEL and the dynamical behavior of glass 
forming liquids has become an active field of research. Among the thermodynamic 
formalisms amenable to numerical investigation of the PEL properties, a central 
role is played by the Inherent Structure (IS) formalism introduced by Stillinger and 
Weber Properties of the potential energy landscape, such as depth, number and 
shape of the basins of the potential energy surface are calculated and used in the 



evaluation of the liquid free energy in supercooled state |18, 19, [2(J, glj]. In the IS 
formalism, the system free energy is expressed as a sum of an entropic contribution, 
accounting for the number of the explored basins, and a vibrational contribution, 
expressing the free energy of the system when constrained in one of the basins JlTj . 

In this manuscript we review the numerical techniques which allow the evaluation 
of the statistical properties of the PEL for atomic and molecular systems [^2[ [01 
^H, HJ, ^l). The paper is organized as follow: Sec. || provides a brief introduction to 



the inherent structure formalism, introduced by Stillinger and Weber 1 17 . Within this 
formalism, an exact expression for the liquid free energy, based on landscape properties 
can be derived. Sec. || reviews the numerical techniques which allow a precise numerical 
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evaluation of the liquid free energy. Sec. [| describes the numerical techniques requested 
for the evaluation of the inherent structure energies. Sec. || discusses techniques 
to evaluate the vibrational component of the free energy. Sec. ^ shows how, from 
the previous information, is possible to evaluate the configurational entropy. Sec. [7] 
discusses a possible modellization of the statistical properties of the landscape, based 
on the hypothesis of a gaussian distribution of basin's depth |2^, |l9[ , and compares 
the predictions of the model with numerical results for a molecular system. 



2. The free energy in the IS formalism 

In the IS formalism , the free energy of a supercooled liquid is expressed in term 
of the statistical properties of the potential energy landscape (PEL). The potential 
energy surface is partitioned in so-called basins, defined as the set of points such that 
a steepest descent path ends in the same local minimum. The configuration at the 
minimum is called inherent structure (IS) and its energy and pressure are usually 
indicated as ejs and Pis- The partition function can be expressed as sum of the 
Boltzmann weight over all the basins, i.e., as sum over the basin partition functions. 
As a result, the Helmholtz liquid free energy F(T, V) can be written as Jl7[ ] : 

F(T,V) = (ei S (V,T)) -TS conf (T,V) + f vib (T,V), (1) 

where 

• (firs) is the average energy of the explored local minima at temperature T and 
volume V; 

• fvib is the vibrational free energy, i.e., the free energy of the system constrained 
in one basin, a quantity depending on the shape of the explored basins; 

• Sconf is the configurational entropy, that counts the number of explored basins. 

The numerical evaluation of F(T, V), (eis(T, V)), and f v a(T, V) is sufficient for 
calculating S CO nf and, from it, the number of basins Q(eis)deis with depth between 
eis and e/s + deis- Indeed, in the thermodynamic limit, lnil(ejs) can be derived 
from a plot of S con f vs (ej<?) (parametric in T). This quantity, together with the (eis) 
dependence of f V ib , provides a precise quantification of the statistical properties of 
the landscape. 



3. Numerical evaluation of F(T, V) 

This section describes the numerical techniques used to evaluate the liquid free energy, 
based on thermodynamic integration [^2], |2(| [18], |2l|]. First, a path in the (T, V") 
plane, connecting the ideal gas state to the desired state point, has to be selected. 
The selected path must avoid the liquid-gas first order line. A convenient choice 
is a constant temperature path (with T = T Q higher than the liquid-gas critical 
temperature) from infinite volume to the desired volume, followed by a constant 
volume path from T a down to the range of temperature of interest. 

In the general case of a system of N rigid molecules, the ideal gas free energy is 



F ig (T,V,N) = -Nk B Ty + -lnTT-lntz + ln 

where A = 27T ™t B , Tt^ = 8,r - with \i denoting x, y, or z 1^ is the inertia 
moment of the molecule with respect the axis fx, and In v accounts for the molecular 



Vy/A 3 K x KyK z 
N 
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symmetry. In the case of C2 V molecules (like water) v = 2, due to the two possible 
degenerate angular orientations of the molecule j27| . 

To perform the thermodynamics integration along the isotherm T , one needs 
to select about 20-30 state points at different volumes (Fig. (a) ) . Of course, the 
smallest chosen volume must coincides with the final volume V . The largest V value 
(Voo) must be chosen in such a way that the vast majority of the molecular interactions 
are binary collisions, i.e., that the volume dependence of the pressure is well described 
by the (first order) virial expansion. At large volumes, although dynamics is very fast, 
care has to be taken to run the simulation long enough to sample a large number of 
binary collisions. 

The free energy at (T Q , V ) can be calculated as 

F(T a , V ) = F lg (T a , V ) - [ ° dVP ex (T , V) + U ^ V °\ (3) 

where U(T Q , V a ) is the potential energy and P ex (T a , V) is the excess pressure, i.e., the 
pressure in excess to the ideal gas pressure. The calculated P ex (T a ,V) curve can be 
fitted according to the polynomial in powers of V^ 1 (Fig. [|-(b)): 

n 

P ex (T ,V) = Y, a ^ V ~ ik+1) > W 
k=l 

giving 

F(T , Vo) = F lg {T , Vo) + £ a ^ T fo~ k + ^L^l . ( 5) 

fe=i ° 
To perform the thermodynamic integration along a constant V path, it is 
necessary to evaluate the internal energy U(T, V ) as a function of T, from T down to 
the lowest state where equilibration of the system is feasible (Fig. [j]-(c)). The resulting 
free energy F(T, V ) can be calculated as 

F(T, Vo) = F(T , Vo) + 3R \n{T/T ) + £ ^L^E^L. (6) 

The 3i?ln(T/T Q ) term accounts for the ideal gas contribution to the free energy. 
Again, a fit of U(T, V ) vs T is requested to evaluate the integral in the above 
expression. One possibility, which has been often found very successful for dense 
systems (small V ) |2g, |2^| , is to fit U(T, V a ) vs T according to the Tarazona law P8| , 
i.e., U(T, V ) = b (V o )+bi(Vo)T 3 / 5 . Of course, for the present purposes, any functional 
form which correctly represents U(T, V a ) can be selected. 

In summary, performing thermodynamic integration, an accurate numerical 
expression for F(V, T) can be obtained. 



4. The average inherent structure energy (ejs) 

This section describes how to calculate the average inherent structure energy 
(eis{T,V)). Recently, it has been shown that, on cooling at constant volume, on 
entering in the supercooled region, the system starts to explore basins of lower and 
lower ej g p2||. The T dependence of the average explored basin depth follows a 1/T 
law |H |19[ [23|, for fragile liquids. Note that for silica, the prototype of a strong 
liquids, the 1/T-law is not observed and (e/s(T, V)) appears to approach a constant 
value on cooling Q. 
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Figure 1. (a) Thermodynamic integration paths used to calculate the total free 
energy at the thermodynamical points of interest, starting from the ideal — non 
interacting- gas state, (b) Excess pressure at T = T as a function of volume. 
The open circles are the MD results. The dashed line is the the first term of the 
virial expansion to the excess pressure; the solid line is a third order polynomial 
fit to the entire set of data, (c) Potential energy (open circles) at the volume V 
in the entire temperature range considered; the solid line is the fit of the data. 
The data are from our simulation |2l| of a system N = 343 molecules modelcd_hy 
the Lewis and Wahnstrom model for Orthoterphenyl |G9|, whose dynamics [ pf| 
and thermodynamics features have been studied in detail. 

In order to evaluate (eis(T, V)) one needs to perform steepest descent potential 
energy minimizations for a statistically representative ensemble of equilibrium 
configurations, to locate their corresponding inherent structure, i.e., local minima. 
For efficiency reasons, the search for the closest local minima is performed using 
the conjugate gradient algorithm [ j34| . In this algorithm, the system evolves along 
a sequence of straight directions until the minimum is reached. In each step, the 
new search direction recalls the directions already explored, improving the algorithm 
efficiency. In rigid molecule systems, each step is composed by a sequence of 
minimization of the center of mass coordinates, followed by a minimization of the 
angular coordinates. Rotations around the principal axis of the molecule are often 
chosen. The minimization procedure is continued until the energy changes less than a 
preselected precision. Since the change in (e/s(T)) in supercooled states is often less 
than one per cent of its own value, a high precision is requested in the minimization 
procedure. 

In Fig. ||we show (eis) (a) as a function of T and (b) as a function of 1/T for a 
rigid molecular model. 
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Figure 2. (erg) as a function of T (a), and as a function of 1/T (b). Data are 
from simulations PjJ of a system TV = 343 molecules modeled by the Lewis and 
Wahnstrom modcHor Orthotcrphcnyl [E9|. 



5. The vibrational free energy 

The vibrational free energy f vlb (e IS , T, V) = U vib (e IS , T, V) - TS vib (e IS , T, V) is the 
free energy associated to the exploration of a basin of depth ejs at temperature T and 
volume V. f V ib{s-iSiT,V) takes into account both the kinetic energy of the system 
and the local structure of the basin with energy eis- From a formal point of view, it 
is defined as 

f^is,T,V) = -k B T\n^- 1 ^^ n(e IS )de IS , 

where ^ is the sum on all the basins with energy depth as- The integration of the 
Boltzmann factor is performed over all points in configuration space associated to the 
selected basin. Here A M = (27r/ Al fc s T) 1 / 2 /h, X = h(2nmk B T)^ 1 ^ is the de Broglie 
wavelength, and I^(r N ) is the potential energy. 

The evaluation of the integral requires the exact knowledge of the shape of the 
PES in the basin and, in general, it will give rise to a complex T-dependence of the 
vibrational energy The best that can be done at the present time is to assume that 
the eis dependence in f V i b is captured by the eis dependence of the density of states 
of the basin, evaluated at the IS configuration |pl|, |3l|]. In other words, the vibrational 
free energy is split into an harmonic contribution (which depends on the curvature 
of the potential energy at the minimum) and an anharmonic contribution, which is 
often assumed basin independent. In molecular systems the Hessian, the matrix of 
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the second derivatives of the potential energy, is calculated numerically selecting as 
molecular coordinates the center of mass and the angles associated with the rotations 
around the three molecular principal inertia axis. Diagonalization is performed with 
standard numerical routines. 

In harmonic approximation, the free energy associated to a single oscillator at 
frequency lo is kBT\n(j3Tiuj). Hence, the basin free energy can be written as 



fvib(eis,T, V) = k B T 



6N-3 

:E 

i=\ 



ln{(3hui(e IS ))y + F anh (T, V), 



with 



and 



U vib (T, V) = (6N-3)^- + U anh (T, V) 



(8) 



(9) 



67V-3 



S mb (e IS ,T,V) = < J2 !- ln 



TiL0 l (e IS ) 



)' + S anh (T,V); (10) 



here the LOi(eis) are the frequencies of the 6N — 3 independent harmonic oscillators 
given by the square root of the 6N — 3 non-zero eigenvalues of the Hessian matrix 
evaluated in the IS. The ()' is the average on all the basins with the same energy e/g. 
Note that the above equations are derived assuming: 



6N-3 6JV _ 

( Hphu}i{e IS )))' = ln(ex P S.=r 



(11) 



For the molecular systems studied so far, this unnecessary approximation introduces 
an error smaller than 1%. The relevant approximation consists in dropping the e/s 
dependence in the anharmonic contribution to the vibrational free energy F an h (and 
of course in U an h(T,V) and S an h{T, V)). In other words, the anharmonicities are 
assumed to be identical in all basins. Under such approximation, U an h{T, V) can be 
calculated from the simulation data as 

k B T 



U anh (T, V) = U(T, V) - (eis(T, V)) - (67V - 3)- 
and it can be well fitted by an expansion in powers of T, starting from T 2 , as 



t/anh(T, V) 



E 

k=2 



c k (V)T k . 



(12) 



(13) 



Correspondingly, S an h (T, V ) can be estimated by thermodynamic integration along 
the isochore between temperatures and T as 

No 



^anh(T,F) 



dT dU(T, V) 



T 



dT 



kc k( v ) T k- 



k=2 



k - 1 



(14) 



By incorporating the anharmonic corrections, which in models of simple fragile liquids 
studied so far are not particularly significant ||^, [H) , a good estimate of the basin free 
energy is obtained. We note on passing that for the cases of network forming liquids, 
anharmonic corrections are relevant p3|, p0| and must be taken into account. 

It is important to note that all the ejs dependence in the average basin free 
energy is carried by the term V(V,eis) = (X)ti~ 3 ^ nLLj i( e is)) ■ A parameterization of 
such quantity as a function of e/s allows to simply connect the basin free energy to 
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Figure 3. V (a) and S con f (b) as a function of eis- 



the basin depth. In all models studied so far |L9, 21, 3^, |3lJ a linear relation between 
basin depth and "basin shape'" V, i.e., 



V = a(V) + b(V)e IS , 
correctly describes the simulation data (Fig. 



(15) 



-(a)). 



6. The statistical properties of the landscape 

In the previous section we have discussed how to relate the basin shape to the basin 
depth. In this section, we exploit the formulation of the liquid free energy in the IS 
formalism to evaluate the number of PEL basins as a function of the basin depth. 

This quantity is of primary interest both for comparing with the recent theoretical 
calculations |jl5| , p6| and to examine some of the proposed relation between dynamics 
and thermodynamics J3(], ^) connecting a purely dynamical quantity like the 
diffusion coefficient to a purely thermodynamical quantity (S con f). The number of 
basins f2 as a function of the basin depth eis has been recently evaluated for a 
few models p8|, p2| , |Ts| , |23| , [l9[ |20| , pi| , pl[ | and from the analysis of experimental 
data |39|, ^llj . S con f {T, V) — the logarithm of the number of explored basins - 
can be calculated as difference of the entropic part of Eq. (||) and Eq. (|To|), i.e., as 

S conf (T, V) = S(T, V) - S mb {T, V). (16) 

In the thermodynamic limit, when fluctuations are negligible, a parametric plot in T of 
Sconf (T, V)/kB vs {eis{T, V)) provides an accurate estimate of the number of basins 
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of depth ejs- This information, together with the information of the e/g dependence 
of the basin shape (or volume) (Eq. (|l5|)) completely defines the statistical properties 
of the landscape, at least in the range of ejs values sampled by the system in the 
studied T region. The availability of S conf (e IS ) (Fig. |-(b)) and V{e IS ) (Fig. |-(a)) 
opens the possibility of a modellization of the thermodynamic of the system in terms 
of landscape properties, as discussed in the next section. 

7. The Random Energy Model: the Gaussian landscape 

A modellization of the statistical properties of the landscape is the next conceptual step 
in the development of a thermodynamic description of the liquid in the IS formalism. A 
possible modellization, which appears to be consistent with the numerical evidence for 
fragile liquids, is based on the hypothesis that the the number fl(ejs)dejs of distinct 
basins of depth between eis and eis + deis m a system of N atoms or molecules is 
described by a Gaussian distribution |^2|, |25|, |23|, |l^, |31j| , i.e., 

-(e IS ~E ) 2 /2<7 2 

n(e IS )de IS = e aN de ls . (17) 

(27T(7 ) / 

Here the amplitude e aN accounts for the total number of basins, E Q has the role 
of energy scale and a 2 measure the width of the distribution. One can understand 
the origin of such distribution invoking the central limit theorem. Indeed, in the 
absence of a diverging correlation length, in the thermodynamic limit, each IS can be 
decomposed in a sum of independent subsystems, each of them characterized by its 
own value of eis- The system IS energy, in this case, will be distributed according to 
Eq. ©. 

Within the two assumptions of Eq. (|j]) — Gaussian distribution of basin depths 
- and Eq. ( p"5| ) — linear dependence of the basin free energy on eis — an exact 
evaluation of the partition function can be carried out. The corresponding Helmholtz 
free energy is given by [jl9| 

F(T, V) = - TS conf (T, V) + (e IS (T, V)) + f mb (E\ T) (18) 

+ k B Tb(V)(< e IS (T,V) > -E'), 

Here E' is an arbitrarily fixed energy, 

(e IS (T)) = (E (V) - b(V)a 2 (V)) - (3a 2 (V), (19) 

and 

S CO nf{T)/k B = a{V)N - ((e IS (T, V)) - E (V)) 2 /2<j 2 {V). (20) 

Note that, from a plot of (eis(T)) vs 1/T, one can immediately evaluate two of 
the parameters of the gaussian distribution, a 2 (from the slope) and E (from the 
intercept). Similarly, from fitting 5 con /(T) according to Eq. (|20j), one can evaluate 
the last parameter a (see Fig. ||). 

The fitting parameters a(V), E Q (V), and cr 2 (V) depend in general on the 
volume. A study of the volume dependence of these parameters, associated to the In- 
dependence of the shape indicators (a and b in Eq. (|l5|)) provides a full characterization 
of the volume dependence of the landscape properties of a model, and offers the 
possibility of developing a full equation of state based on statistical properties of the 
landscape. 
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When comparing numerical simulation data and theoretical predictions, - 
Eqs. ( |l9l ) and ( po| ) the range of temperatures must be chosen with great care. Indeed, 
at high T, the harmonic approximation will overestimate the volume in configuration 
space associated to an inherent structure. While in harmonic approximation such a 
quantity is unbounded, the real basin volume is not. Indeed, the sum of all basins 
volumes is equal to the volume of the system in configuration space. Anharmonic 
corrections, if properly handled, should compensate such overestimate, but at the 
present time, no model has been developed to correctly describe the high T limit of 
the anharmonic component. Numerical studies have shown that the range of validity 
of the present estimates of the anharmonic correction does not extend beyond the 
temperatures at which the system shows already a clear two-step relaxation behavior 
in the dynamics. Indeed, the presence of a two-step relaxation is a signature that the 
system spends a time larger than the microscopic characteristic times around a well 
defined local minima. 

8. Conclusions 

In this paper we have discussed the numerical techniques employed to evaluate the 
statistical properties of the PEL for molecular systems. These numerical calculations 
are limited to the region of temperatures and volumes where equilibrium configurations 
can be numerically generated. Still, very simple arguments can be presented which 
allow to generalize the results and formulate a full thermodynamic description of the 
supercooled liquid state, only in terms of the statistical properties of the PEL. 

The possibility of partitioning the free energy and its thermodynamic derivatives 
as a sum of configurational and vibrational degrees of freedom has been recently 
exploited to derive a satisfactory description of the equation of state ^| |3l| for 
supercooled liquids only in terms of PEL properties. A better understanding of 
the nature of each contribution (configurational and vibrational) to quantities as 
the total pressure of the system is achieved. At the same time, the availability of 
detailed estimates for the landscape properties strongly suggests a generalization of 
this approach to out of equilibrium conditions. It has been recently shown [[iif that 
if the system ages exploring the same basins visited in equilibrium, it is possible to 
write an out of equilibrium equation of state expressing P as a function not only of V 
and T but also as a function of the (time dependent) depth of the explored basin. 

The availability of numerical estimates for the statistical properties of the PEL 
in models of simple liquids should encourage theoreticians to develop schemes for the 
analytic evaluation of these quantities. If this goal were reached, the understanding 
the thermodynamics of supercooled liquids and glasses would improve significantly. 
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